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The equations in conservative form for nonlinear waves modeling on a liquid film flowing down 
a vertical plane have been investigated. It has been found that in the computational domain 
extended along the transverse axis the equations with boundary conditions are invariant under 
parity transformation. It is numerically shown that for moderate Reynolds numbers the steady- 
state travelling solutions of the equations have the detected symmetry. It is demonstrated that 
using this symmetry for the numerical solution of the problem by Galerkin methods significantly 
increases the efficiency of calculations. 
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I. INTRODUCTION 


Film flows are widely used in heat and mass transfer 
technologies. Theoretical study of the problem started 
in 1916 with Nusselts work that has become fundamen¬ 
tal for the theory of film condensation. The research 
provided the exact solution for waveless flowing of a thin 
layer of viscous liquid over a solid inclined wall. The nat¬ 
ural development of waves on the free surface of a film 
additionally intensifies the process of heat transfer from 
the heated substrate. The study of the basic characteris¬ 
tics and the form of these waves in the pioneer works of 
Kapitza [li] laid the foundation for a new scientific direc¬ 
tion. 

Many of the researchers (for example, i: i) use low¬ 
dimensional Galerkin methods to derive reduced systems 
of equations for coefficients at respective basis functions, 
depending on the transverse coordinate. The main ad¬ 
vantage of such approaches results from rather small 
number of equations (two or three) at a lucky choice 
of the basis. In our work we point to a specific sym¬ 
metry of governing equations and a wide class of their 
solutions, including those of interest from experimental 
point of view. It allows two-fold reduction of the num¬ 
ber of equations and gives a deep understanding of suc¬ 
cessfulness of the known models, derived earlier by ex¬ 
trapolation of analytical results, obtained only for small 
Reynolds numbers. 

The full problem formulation for waves on isothermal 
flowing film is clear and includes a system of Navier- 
Stokes and continuity equations with appropriate bound¬ 
ary conditions at the wall and at the free surface. This 
formulation implies the solution of the problem in a 
changing flow area unknown in advance, which greatly 
complicates the mathematical and numerical simulation. 
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One of the approaches to solving the moving boundary 
problem appeared in the mid 80-ies of the last century. 
It comes to rewriting of hydrodynamic equations in new 
variables, transforming the flow area into the strip of 
constant thickness (Fig. [T|: 

x = x, r\= , V , t = t (1) 

n(x, t) 

Here h is an instant local film thickness. The coordinate 
system m is non-orthogonal, so the normal vector for¬ 
mulation of equations is inapplicable. For this reason, 
many authors (see for example fli) reduce the method 
to a simple change of variables without transformation of 
vectors and tensors, contained in the original equations. 



FIG. 1. The flow area in x-y and x-g coordinates 


For the case of a free falling film on a vertical plane, the 
transformation m was performed in [(5] for the hydrody¬ 
namic equations written in tensor form invariant under 
arbitrary coordinate transformation. As a result, in the 
long wavelength approximation the following system was 
obtained: 
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dh d(hu) d(hv) 

dt dx dr] 


(2b) 


Here a is a surface tension, p is a density, v is a kine¬ 
matic viscosity, g is a free fall acceleration, u and v are 
contravariant components of velocity corresponding to 
the coordinates x and 77, respectively. It is clear that 
since x has remained an original Cartesian coordinate, 
the component u(x,r],t) in contrast to v(x,r],t) remains 
the longitudinal component of physical velocity. 

Boundary conditions on the solid wall are: 


77(3;, 0 , f) = 0 , v(x, 0 ,t )=0 ( 3 ) 

No-stress Cal) and kinematic conditions ( 14 bl) on the free 
surface are: 


v(x , ??', t) —>• —v(x, —77', t) ( 5 c) 

It means in particular that there are two types of solu¬ 
tions of these equations in the extended strip 77' £ [—1,1], 
Solutions of the first type are characterized by symmetry: 
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It is clear that the solutions of the first type satisfying 
the no-slip and no-penetration boundary conditions on 
both boundaries: 

u(x, — l,t) = v(x, — l,t) = 0 ( 7 a) 


|^(a 7 , 1 ,t) = 0 

( 4 a) 

77(37,1 , t) = 0 

( 4 b) 


In this problem there are three unknown functions, one 
of which h is the function of only two variables (x, t). In 
the calculated domain ( x , 77) there are 2 unknown func¬ 
tions - u(x,r],t) and 77(37, 77, f), therefore, 2 equations are 
necessary to find them. In the long-wave model these two 
equations are: continuity equation and momentum bal¬ 
ance in the longitudinal direction. A kinematic condition 
on the free surface (77 = 1), allowing finding the function 
h(x, t). In contrast to usual presentation of the kinematic 
condition on a free boundary in Cartesian coordinates, 
where it is represented as a differential equation: 

/ , . dh dh 

v{x,y = h,t) = — + u{x,y = h,t)— 

in curvilinear coordinates ( 1 ) it takes the form ( 14 bl) : 


II. SYMMETRY 


u(x, 1 , t) = v{x, 1 , t) = 0 ( 7 b) 


are the solutions to original problem ©-© in a half-strip 
[— 1 , 0 ]. Indeed on the boundary 77' = —1 the conditions 
© are performed, and the boundary condition of no¬ 
penetration at 77' = 0 (kinematic condition v(x , 0, t) = 0) 
is satisfied automatically, due to fact that the contravari¬ 
ant transverse component of the velocity v for the first 
type solutions is an odd function with respect to 77'. And 
since the function u is even, the dynamic condition on 
this boundary is automatically performed as well: 


du 

dip 


7(^,0 ,t) = 0 


Solutions of the second type do not have such symme¬ 
try, but in virtue of ©, if there is a solution ui(x,r]',i) 
and v\(x, 77', t), then there is a solution: 


u 2 (x,r)',t) = ui(x, -77', t) 


At that 


772(37,77', f) = —77i (x, -77', t) 


Ul(x, 77', t) / 77 i(x, -77', t) 


Let new transverse coordinate be 
t/ = ?7 - 1 

so that the flow area lies in the interval 77' £ [—1,0]. 
Notice that 77' is expressed in the original Cartesian co¬ 
ordinates as follows: 

, = y _ 

^ h(x,t) 

It is easy to see that the equations © are invariant under 
the transformation: 


/ . / 

T) -T) 

( 5 a) 

73(3;, 77', t) —>• u(x, —77', t) 

( 5 b) 


Vi(x,T]',t) ± — 77 !(37,-77', f) 

However, if the second type solutions for the problem 
©-© are extended into interval 77' £ [0,1] then, bound¬ 
ary conditions m at 77' = 1 are not to be performed. 

Considering higher orders of smallness makes the prob¬ 
lem as a whole and boundary conditions on free surface 
© in particular much more complicated. All this will 
most probably lead to symmetry © breaking. Anyway, 
the pointed-out symmetry will be destroyed even at the 
used approximation if other conditions at free boundary, 
for example, the ones taking into account gas influence 
on film flowing, are considered. Nevertheless, the found 
symmetries for the model of a boundary layer, are quite 
interesting themselves. 

The above symmetry properties of equations may be 
useful for finding solutions to the problem ©-©. 
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III. CALCULATIONS 


Hence, assuming a n = 1, for e obtain: 


For numerical solution the problem ©-0 was written 
in a dimensionless form: 


e = 


6 3Re 5/3 
53VW/3 


(9) 


eRe 


dQ d (Q^\ d ( 

dt dx\h) drf \ 
1 d 2 Q 3 1 / 3 Fi V 3 

h 2 dr\' 2 Re 5 / 3 



e 3 Reft 


d 3 h 

fa 3 


(8a) 


dh dQ dV_ 
dt dx dr]' 


(8b) 


Q(x,—\,t) = V{x,—l,t)=Q (8c) 


^(s,0,t) = 0, V(x,0,t)=0 (8d) 

We introduce new functions Q = uh , V = vh relative to 
which the equation (I2bl) is linear. Here Re is Reynolds 
number, Fi = a 3 /p 3 gv^! 3 is a film number characterizing 
physical properties of liquid, £ is a long-wave parameter, 
Mo, lo, ho and Io/uq are characteristic scales of velocity, 
length, thickness and time, respectively. The solution of 
the linear problem ©-0 shows that the unperturbed 
flow 


h 0 = 1, Qo(v') = |(1 - rf 2 ), V 0 = 0 

is unstable with respect to linear perturbations of the 
type 

[Q\ h', V'} = [<X, 1, U a '] {rf) h! a exp (ia(x - ct )) + c.c. 

for an interval of wave numbers a less than a n . Here 
c = c r + ici is complex phase velocity. The scales of 
length and thickness (and hence the parameter e) may be 
selected so that the wave number of neutral perturbations 
(ci = 0) for any values of Reynolds numbers was equal 
to a n = 1. With this normalization for the undisturbed 
flow the area of instability (c* > 0) with respect to linear 
perturbations would take an interval 0 < a < a n = 1. 
Accordingly, the wave numbers of stable perturbations 
would be in the area a > a n = 1. For such normaliza¬ 
tion it is necessary to determine the relationship between 
parameters e, Re, Fi from the solution to a linear prob¬ 
lem. For the problem ©-0 such relationship in the 
general case cannot be determined analytically. 

In fj] it is shown that at low Reynolds numbers, the 
problem ©-0 are reduced to the equation of Kuramoto- 
Sivashinsky (Nepomniashchi), where the wave number of 
neutral perturbations is known (see for example [§j: 

2 2 6 3Re 5/3 

£ ~ a " ~ 531 / 3 ^ 1/3 


This is the ratio that was used to consider the problem © 
in the area of moderate Reynolds numbers. In accordance 
with calculations, up to values eRe = 6 (for water films 
is Re ~ 20) a n differs from unity by less than 0.002. 
Taking into account expression ©, the coefficient of the 
capillary term in the equation (I8al) takes the form: 


3 i/ 3 F i /3 

Re 5 / 3 


e 3 Re 



so only parameter eRe remains in the problem. To solve 
a linear problem, the algorithm described in j§| was used: 
the equations © was reduced to one equation of the mod¬ 
ified stream function, then the equation was linearized 
and solved by the shooting method from the substrate 
rf = — 1 (to solve the Cauchy problem we used the 4th- 
order Runge-Kutta method). The dynamic condition on 
the free boundary determined a free coefficient in a homo¬ 
geneous solution to the equation. The phase velocity was 
found by Newton’s method from the kinematic condition 
on the surface. According to direct numerical calcula¬ 
tions, extending the linear solution to the problem © 
into the interval rf G [0,1] we obtain solutions with sym¬ 
metry both for stable and unstable perturbations. Two 
examples of such calculations are shown in Fig. [2] 
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FIG. 2. Real (Q r ,V r ) and imaginary ( Qim,Vim) parts of Q 
and V depending on dimensionless transverse coordinate in 
extended flow area. 

The growth of unstable linear perturbations may be 
eventually stopped by the action of nonlinear effects, 
resulting in the formation of steady-state non-linear 
regimes. Therefore, in the studies of wave flows of liquid 
films much attention is paid to the steady-state travel¬ 
ling waves. In this work we will limit consideration to 
periodic steady-state travelling solutions of system ©: 

[Q, h, V ] = [Q, h, V] (x, rf ), x = a(x - ct) 
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In connection with the results obtained for linear solu¬ 
tions, the following question arises, how wide is the class 
of steady-state travelling solutions of the problem ©-©, 
for which the symmetry © is performed as well? 

Currently the interest in Galerkin spectral methods as 
applied to the problems of fluid motion has markedly in¬ 
creased. To describe film flows a number of new model 
systems obtained by this method were introduced. The 
minimum number of members in Galerkins approaches, 
required for a good approximation of solutions, greatly 
depends on how well the basis functions were chosen. For 
the numerical study we used the pseudospectral method, 
in which functions depending on the transverse coordi¬ 
nate rf were expanded in series on Chebyshev polynomi¬ 
als Tf. 


Q = Qi( x )Ti{r]') 

( 10 a) 

V = Y,V i {x)T i (r 1 ') 

( 10 b) 

and functions depending on the longitudinal coordinate 
were expanded in Fourier series: 

Q I (i) = ^Qfe ifc£ 

( 11 a) 

Vi(x) = J2 v f eik * 

(lib) 

h(x) = Hkeiki 

( 11 c) 


Limiting to the first N harmonics in (fTTh and substi¬ 
tuting expressions (fl 0 l) - (fTTl) in system © written for a 
variety of Chebyshev nodes: 

^ = C0S (^Ar^ J = 1 - M 

and satisfying the boundary conditions at the points 
rf = —1 and rf = 0 , we come to a finite-dimensional 
system of nonlinear algebraic equations for coefficients 
V k 1 H k ,c. The resulting problem was solved by the 
iterative method of Newton-Kantorovich: 

[Q, h, V , c, a] n+1 = [Q, h, V, c, a] n + [8Q, Sh , 8V, 8c , Sa] 

Here the symbol 8 indicates the variation of the corre¬ 
sponding value. After linearization relative to variations 
the problem is reduced to a system of linear algebraic 
equations. 

Using algorithm described above a series of calcula¬ 
tions of steady-state travelling solutions of the system 
© in the strip 77 ' £ [—1,1] was performed. The param¬ 
eter eRe varied in the range [0.03 8 ]. For the case of 

water this means a change of Reynolds number in the 
range [1 4- 30]. A substantial part of the calculations was 
performed on the bases msD with the number of polyno¬ 
mials Ti - 9, 11, 13. The calculations have shown that 


the difference in the waveform obtained on the bases m, 
containing 9 and 11 polynomials, was discernible in ar¬ 
eas with large gradients of thickness (on the wave front 
and in the field of precursor). But even in these areas 
the results, built on the bases of 11 and 13 polynomials, 
coincided with graphic accuracy. The most remarkable 
result is that the values of coefficients obtained during 
calculations for all odd Chebyshev polynomials turned 
out to be equal to null with machine precision. Thus, it 
has become clear that all found steady-state travelling so¬ 
lutions of system ©-© are characterized by symmetry 
©■ 

Considering these results, for comparison we have ap¬ 
plied an alternative method of solving the system ©, 
where for the function Q we used only the even Cheby¬ 
shev polynomials, and for the function V - only the odd 
ones: 

q = (*) c 7 * (»?') -!) ( 12a ) 

v = ^i+1 ( 5 ) (T«+1 {rf) - rf) (12b) 

Such a choice provides automatic performance of bound¬ 
ary conditions at the interface and on the solid wall. Of 
course now it is sufficient to consider points rf only in a 
half-band [—1,0]. This representation gives the solution 
to our problem ©-([3]). Several control calculations car¬ 
ried out in the basis ED. consisting of 3 even Chebyshev 
polynomials for the longitudinal velocity, have shown ex¬ 
cellent agreement with the results, obtained on the full 
basis ED from 9 polynomials. 

The obtained steady-state travelling solutions of the 
problem ©-(© are adequately presented in Here we 
give one of these solutions as an example. Fig. [3] presents 
the form of the surface for this case. For a more detailed 
description of the solution, it is useful to present profiles 
of functions Q and V in different wave sections. The 
most interesting are areas with large thickness gradients 
(leading wave front and precursor). Here, for example, 
the profile Q is significantly different from the parabolic, 
and it may contain a point of inflection. In these areas 
reverse flows may occur as well. Sections in enlarged part 
of Fig. [3] suitable from this point of view are numbered 
with points A-E. The corresponding profiles of functions 
Q and V are shown in Fig. []] and Fig. [5] To demon¬ 
strate symmetries © profiles are given in the interval 
rf £ [-1,1]. 

Of course, on the basis of performed calculations it 
cannot be considered strictly proved that all possible 
steady-state travelling solutions of the problem ©-© 
in the band rf G [—1,1] have symmetry ©, but such a 
conclusion is very likely. At least for a special case of 
periodic steady-travelling modes having positive solitons 
as limits, which is important from the point of view of 
the experiment, we can assume that this is true. There¬ 
fore, we can conclude about the usability of shorter and 
convenient basis functions, based on symmetry © of the 
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FIG. 3. Profile of the wave 



FIG. 4. Profiles of the function Q for selected wave sections 
A - E 


equations ©, for calculating the steady-state travelling 
modes of wave motion of the film. 

For example, the symmetry is specific for the function, 
used in the well-known model of Shkadov. In this context 
let us mention the model, used in the well known series of 



FIG. 5. Profiles of the function V for selected wave sections 
A - E 


articles by Ruyer-Quil and Manneville and show 

that there the basic functions are also symmetric. To 
represent the longitudinal velocity u the authors used a 
basis consisting of three polynomials: 

1 2 

9o = V~ 2 11 


17 2 7 3 7 4 

9i=V~ Y V + —T) 


13 , 
92 = 9- y'7 


57 


111 


„ 99 5 

-r) H- rf 

' 16 ' 


33 

32 ? 


The choice of such a basis was almost pushed by the an¬ 
alytical results, obtained in the case of weakly-nonlinear 
perturbations at small flow rates, when the Reynolds 
numbers are of the order of unity. In this basis the ve¬ 
locity is written as follows: 


2 

u(x, 77 , t) = b o( x ’ t )9j(v) 

0=0 


By means of shift transforming g —»• if the basis g 3 of 
Ruyer-Quil, Manneville takes the following form: 


__ 1 1 
90 2 2 9 


1 , 2 /2 7 /4 

91 = -12 + 3 V ~ 12 ?? 


- 1 _ 1£ '2 , 51 /4 _ 33 /6 
92 32 32^ + 32 V 32 V 


Now it has become clear that the polynomials g 3 con¬ 
tain only even terms relative to the coordinate if. There¬ 
fore, all possible expressions of the longitudinal velocity 
presented with this set will be symmetric with respect to 
the coordinate rf. 
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IV. CONCLUSION 

The study equations in conservative form for simulat¬ 
ing the dynamics of nonlinear waves on the liquid film 
flowing down a vertical plane has been performed. It has 
been found that in the computational domain r\' £ [—1,1] 
extended along the transverse coordinate the equations 
© with boundary conditions ([l])-© is invariant under 


parity transformations ©. It was numerically proved 
that at moderate Reynolds numbers the steady-state 
travelling solutions of the problem (f2 |l -([4 jl in the area, 
extended along the transverse coordinate, are character¬ 
ized by the detected symmetry ©• Relevance of this 
symmetry consideration for the numerical solution of the 
problem using spectral and pseudospectral methods has 
been demonstrated. 
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